for (package in c("tidyverse","fpp3", "GGally", "normtest")) {
if (!require(package, character.only=T, quietly=T)) {
install.packages(package)
library(package, character.only=T)
}
}
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.5 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.0.2 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## ── Attaching packages ──────────────────────────────────────────── fpp3 0.4.0 ──
## ✓ lubridate 1.8.0 ✓ feasts 0.2.2
## ✓ tsibble 1.1.1 ✓ fable 0.3.1
## ✓ tsibbledata 0.4.0
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## x lubridate::date() masks base::date()
## x dplyr::filter() masks stats::filter()
## x tsibble::intersect() masks base::intersect()
## x tsibble::interval() masks lubridate::interval()
## x dplyr::lag() masks stats::lag()
## x tsibble::setdiff() masks base::setdiff()
## x tsibble::union() masks base::union()
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
# Contraste para los coeficientes
my_t_test <- function (object, ...)
{
par <- rbind(t_stat=tidy(object)$statistic, p_value=tidy(object)$p.value)
colnames(par) <- tidy(object)$term
if (NCOL(par) > 0) {
cat("\nt-test:\n")
coef <- round(par, digits = 4)
print.default(coef, print.gap = 2)
}
}
# Gráfico de correlogramas de residuos
my_tsresiduals <- function (data, ...) {
if (!fabletools::is_mable(data)) {
abort("gg_tsresiduals() must be used with a mable containing only one model.")
}
data <- stats::residuals(data)
if (n_keys(data) > 1) {
abort("gg_tsresiduals() must be used with a mable containing only one model.")
}
gg_tsdisplay(data, !!sym(".resid"), plot_type = "partial",
...)
}
# Validación cruzada anidada
nested_cv <- function(df, h, last_train, string_formula){
nested_errors <- vector()
for (i in seq(last_train, last(df$date), h)){
train <- df %>% filter(date<=i)
test <- df %>% filter(date>i)
fitted_model <- train %>%
model(arima = ARIMA(as.formula(string_formula)))
h_forecast = min(dim(test)[1], h)
fc <- fitted_model %>%
forecast(h=h_forecast)
test_err <- fc %>%
accuracy(test) %>%
select(MAPE)
nested_errors <- c(nested_errors, test_err$MAPE)
NewList <- list("errors"=nested_errors, "mean"=mean(nested_errors))
}
return(NewList)}
# Test de autocorrelacines de los residuos
autocorrelation_test_plot <- function(aug, dof = 4, m = 7, h = 5, alpha = 0.05){
vec <- c()
num_lags = seq(1, h*m)
for (i in num_lags){
vec <- c(vec,aug %>% features(.resid, ljung_box, lag=i, dof=dof) %>% .$lb_pvalue)
}
autocorr_pvalues_resid <- tibble(
lag = num_lags,
p_value = vec,
incorelated = p_value >= alpha
)
plot <- autocorr_pvalues_resid %>%
drop_na() %>%
ggplot(aes(lag, p_value, color = incorelated)) +
geom_point() +
geom_hline(aes(yintercept = alpha), linetype="dashed", color = "indianred2")
newList <- list("values" = autocorr_pvalues_resid, "plot" = plot)
return(newList)
}
Antes de empezar a hacer el esco2dio hay que distinguir la parte de los datos que se utilizará para construir la fórmula (para modelizar) de aquella que se utilizará para validar los resultados y cuyos datos, no deben ser utilizados.
Convertimos la serie en un objeto tsibble, para trabajar de manera cómoda. Analizamos visualmente la serie y las marcas de tiempo para saber cómo dividir los datos.
co2 <- read_csv('co2_mm_mlo.csv')
## Rows: 746 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (3): year, month, value
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
co2$date <- with(co2, sprintf("%d-%02d", year, month)) # juntar mes y año
co2$date = tsibble::yearmonth(co2$date)
co2 = as_tsibble(co2, index = date)
co2 %>%
autoplot(value) +
labs(title = "CO2") +
xlab("Fecha") + ylab(" ") +
geom_point(aes(y = value), size = 1.2, shape = 16, color = "black")
Dado que tenemos datos mensuales y visualmente se aprecia una estacionalidad anual, reservamos el último año para test y el resto para train.
co2_train <- co2 %>% filter_index(. ~ "2005-12-01")
co2_test <- co2 %>% filter_index("2006-01-01" ~ .)
Se comienza la fase de identificación donde se deciden las transformaciones a realizar y el propio ajuste de la serie. Aunque veremos que cuando se realiza la selección de los parámetros del modelo ARIMA, la iteración con la estimación y el contraste son constantes hasta converger a un modelo válido.
Antes de analizar la serie mediante test estadísticos, analicemos gráficamente más en detalle, para confirmar la estacionalidad detectada.
co2_train %>%
gg_season(value, period = "year",labels = "left")
Estacionalidad anual: la forma de las series anuales son iguales.
Descomposición STL para analizar la existencia de tendencia y estacionalidad.
co2_train %>%
model(seats = feasts:::STL(value)) %>%
components() %>%
autoplot()
Se confirman visualmente la estacionalidad anual y una tendencia creciente y luego decreciente.
A continuación, esco2diamos si la serie es estacionaria en media y varianza y acco2amos en consecuencia.
Analicemos qué transformación Box-Cox puede aplicarse para estabilizar la varianza.
lambda <- co2_train %>%
features(value, features = guerrero) %>% pull(lambda_guerrero)
lambda
## [1] 0.7022992
co2_train %>% autoplot(box_cox(value,lambda)) +
labs(y = "Box-Cox transformed TU")
Dado que la serie visualmente muestra una varianza más o menos constante, que el valor Box Cox es más próximo a 1 que a otros valores de transformaciones usuales (logaritmo, raíz cuadrada) y que no ha habido un cambio muy llamativo con la transformación, no realizaremos ningún cambio a priori. Si en el proceso de estimación no conseguimos un modelo adecuado, volveríamos a este punto para valorar de nuevo la transformación.
La serie muestra una tendencia creciente y luego decreciente, es decir, la media va variando a lo largo de la serie. Por tanto, es esperable que se necesite realizar al menos una diferencia para estabilizar la media. Formalmente, se comprueba mediante una prueba de raíces unitarias.
co2_train %>%
features(value, unitroot_kpss)
El p-valor es menor que 0.05, lo que indica que la hipótesis nula es rechazada. Es decir, hay evidencias estadísticas para rechazar la hipótesis nula de que la serie sea estacionaria. Por tanto, habría que realizar una diferencia regular.
co2_train %>%
features(difference(value, 1), unitroot_kpss)
Una vez realizada la diferencia regular, el p-valor del test indica que no hay evidencias estadísticas para rechazar la hipótesis de que la serie sea estacionaria en media.
Graficamos la serie diferenciada y ya observamos que el nivel de la serie no cambia:
co2_train %>% autoplot(difference(value, 1)) + labs(title = "Serie con una diferencia regular")
Debido a la estacionalidad de los datos, quizá deba realizarse una diferencia estacional. Para evaluarlo podemos recurrir a la función unitroot_nsdiffs(), esta evalúa el estadístico fuerza estacional \(F_S\) y sugiere una diferencia estacional si esta es mayor que 0.64.
co2_train %>%
mutate(turnover = difference(value,1)) %>%
features(turnover, unitroot_nsdiffs)
El método sugiere una diferencia estacional. En cualquier caso, para asegurar que la serie es estacionaria, en la fase de identificación del modelo debemos asegurarnos de que las raíces de la parte autoregresiva están fuera del círculo unitario (y, por tanto, 1/raíz dentro del círculo).
En los pasos previos se ha sugerido una diferencia regular y una estacional. Empezaremos teniendo en cuenta con la diferencia regular. Es decir, nos dirigimos a identificar los órdenes \(p, q, P, Q\) de un SARIMA (p,d,q)x(P,D,Q)\(_{12}\) donde ya sabemos que \(d=1\). Vemos que los residuos contienen mucha información.
fit <- co2_train %>%
model(arima = ARIMA(value ~ 1+pdq(0,0,0) + PDQ(0,0,0, period = 12)))
fit %>% my_tsresiduals(lag_max = 36)
Claramente \(p=1\).
Comencemos ajustando la parte estacional mediante un SARIMA (1,0,0)x(0,0,0)\(_{12}\).
fit2 <- co2_train %>%
model(arima = ARIMA(value ~ 1+pdq(1,0,0) + PDQ(0,0,0, period=12)))
fit2 %>% my_tsresiduals(lag_max =36)
report(fit2)
## Series: value
## Model: ARIMA(1,0,0) w/ mean
##
## Coefficients:
## ar1 constant
## 0.9991 0.3188
## s.e. 0.0012 0.0232
##
## sigma^2 estimated as 1.479: log likelihood=-893.43
## AIC=1792.86 AICc=1792.91 BIC=1805.8
my_t_test(fit2)
##
## t-test:
## ar1 constant
## t_stat 853.9515 13.7145
## p_value 0.0000 0.0000
Este valor indica que hay que diferenciar (como ya habíamos comprobado).
Diferenciamos, \(d=1\). SARIMA (0,1,0)x(0,0,0)\(_{12}\).
fit3 <- co2_train %>%
model(arima = ARIMA(value ~ 1+pdq(0,1,0) + PDQ(0,0,0, period=12)))
fit3 %>% my_tsresiduals(lag_max =36)
report(fit3)
## Series: value
## Model: ARIMA(0,1,0) w/ drift
##
## Coefficients:
## constant
## 0.1159
## s.e. 0.0515
##
## sigma^2 estimated as 1.463: log likelihood=-886.2
## AIC=1776.4 AICc=1776.43 BIC=1785.03
my_t_test(fit3)
##
## t-test:
## constant
## t_stat 2.2518
## p_value 0.0247
La ACF muestra una estructura AR estacional \(P=1\). Es decir, SARIMA (0,1,0)x(1,0,0)\(_{12}\).
fit4 <- co2_train %>%
model(arima = ARIMA(value ~ 1+pdq(0,1,0) + PDQ(1,0,0, period=12)))
fit4 %>% my_tsresiduals(lag_max =36)
report(fit4)
## Series: value
## Model: ARIMA(0,1,0)(1,0,0)[12] w/ drift
##
## Coefficients:
## sar1 constant
## 0.9436 0.0068
## s.e. 0.0124 0.0131
##
## sigma^2 estimated as 0.1652: log likelihood=-297.76
## AIC=601.52 AICc=601.56 BIC=614.45
my_t_test(fit4)
##
## t-test:
## sar1 constant
## t_stat 76.1768 0.5210
## p_value 0.0000 0.6026
Introducimos \(Q=1\) (se salen los múltiplos de 12 en la PACF y la correlación más alta en la ACF es la de 12). Es necesario también hacer una diferencia estacional.
SARIMA (0,1,0)x(0,1,1)\(_{12}\).
fit5 <- co2_train %>%
model(arima = ARIMA(value ~ 1+pdq(0,1,0) + PDQ(0,1,1, period=12)))
fit5 %>% my_tsresiduals(lag_max =36)
report(fit5)
## Series: value
## Model: ARIMA(0,1,0)(0,1,1)[12] w/ poly
##
## Coefficients:
## sma1 constant
## -0.9200 0.0024
## s.e. 0.0204 0.0015
##
## sigma^2 estimated as 0.09687: log likelihood=-139.1
## AIC=284.2 AICc=284.24 BIC=297.06
my_t_test(fit5)
##
## t-test:
## sma1 constant
## t_stat -45.108 1.6117
## p_value 0.000 0.1076
Habría que meter \(p=1\), \(q=1\). SARIMA (0,1,1)x(0,1,1)\(_{12}\).
fit6 <- co2_train %>%
model(arima = ARIMA(value ~ 1+pdq(0,1,1) + PDQ(0,1,1, period=12)))
fit6 %>% my_tsresiduals(lag_max =36)
report(fit6)
## Series: value
## Model: ARIMA(0,1,1)(0,1,1)[12] w/ poly
##
## Coefficients:
## ma1 sma1 constant
## -0.346 -0.8874 0.0024
## s.e. 0.044 0.0216 0.0012
##
## sigma^2 estimated as 0.08865: log likelihood=-112.16
## AIC=232.33 AICc=232.4 BIC=249.49
my_t_test(fit6)
##
## t-test:
## ma1 sma1 constant
## t_stat -7.8643 -41.1472 2.0855
## p_value 0.0000 0.0000 0.0375
gg_arma(fit6)
Para completar la fase de contraste evaluamos los residuos. Veamos primero el histograma.
aug <-fit6 %>% augment()
# Histogram
aug %>%
ggplot(aes(x = .resid)) +
geom_histogram(bins = 50) +
ggtitle("Histogram of residuals")
Contraste t-student para contrastar si la media es 0.
# Student's t-Test for mean=0
t.test(aug$.resid)
##
## One Sample t-test
##
## data: aug$.resid
## t = 0.0027261, df = 551, p-value = 0.9978
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.02451750 0.02458564
## sample estimates:
## mean of x
## 3.407357e-05
Como p-value > 0.05 no existen evidencias significativas para rechazar la hipótesis nula de que la muestra tenga media 0.
A continuación comprobamos que los residuos están incorrelados.
# Ljung-Box autocorrelation
aug %>% features(.resid, ljung_box, lag=12, dof=3)
Podemos concluir que los residuos están incorrelados.
Para contrastar la heterocedasticidad se puede utilizar una regresión media-dispersión. Calculando los grupos de manera anual, dado que esa es la estacionalidad de la serie.
log_log <- aug %>% as_tibble() %>%
group_by(year(date)) %>%
summarize(mean_resid = log(mean(.resid+1)), std_resid = log(sd(.resid+1)))
summary(lm(std_resid~mean_resid, log_log))
##
## Call:
## lm(formula = std_resid ~ mean_resid, data = log_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.25683 -0.14625 0.01988 0.20569 0.39826
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.27734 0.04477 -28.530 <2e-16 ***
## mean_resid -0.78133 0.71330 -1.095 0.279
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3035 on 44 degrees of freedom
## Multiple R-squared: 0.02655, Adjusted R-squared: 0.004422
## F-statistic: 1.2 on 1 and 44 DF, p-value: 0.2793
El p-valor de log(media) es mayor que 0.05, no hay evidencias para rechazar la hipótesis de que sean homocedásticos.
Para contrastar la normalidad realizamos el test Jarque-Bera para evaluar la normalidad.
# Jarque Bera test
jb.norm.test(na.omit(aug$.resid))
##
## Jarque-Bera test for normality
##
## data: na.omit(aug$.resid)
## JB = 9.0477, p-value = 0.016
Como el p-valor <0.05 se rechaza la hipótesis nula de normalidad de los residuos.
Finalmente evaluamos la capacidad predictiva del modelo. Calculamos las predicciones en el conjunto de test, calculamos sus errores y los comparamos con los residuos.
# Residual accuracy
resids <- fit6 %>%
accuracy() %>%
select(-c(.model, .type, ME, MPE, ACF1, RMSSE)) %>%
mutate(Evaluation='Training')
# Forecasting
fc <- fit6 %>%
forecast(h=12)
test_err <- fc %>%
accuracy(co2) %>%
select(-c(.model, .type, ME, MPE, ACF1, RMSSE)) %>%
mutate(Evaluation='Test')
# Show errors together
bind_rows(test_err, resids) %>% select(Evaluation, everything())
Además de evaluar los errores con estas medidas globales, podemos evaluar los errores cometidos año a año visualmente. Así, podemos detectar outliers o efectos de calendario no detectados. Pudiendo así corregirlos con modelos más complejos mediante variables exógenas.
aug %>% ggplot() +
geom_line(aes(x = date, y = .fitted), color="red") +
geom_line(aes(x = date, y = value), color="gray24") +
ggtitle("SARIMA train fitted values") +
xlab('Dates') +
ylab(' ') #+ facet_wrap(vars(year(date)), scales = 'free')
También podemos dibujar las predicciones con los intervalos de confianza (aunque en este caso no sean representativos, dada la no normalidad de los residuos). Evaluamos así también la capacidad de generalización del modelo, comparando las predicciones con el conjunto de test.
h <- dim(co2_test)[1]
demanda_plot <- co2 %>% filter(date>last(co2_train$date)-14 & date< last(co2_train$date) + 120 )
fit6 %>%
forecast(h=120) %>%
autoplot(demanda_plot)
Como sabemos que los modelos SARIMA no predicen bien más alla del orden más alto disponible, en este caso 12, los errores en el conjunto de test más allá de 12 datos se dipararían.